Exploring Defective Eigenvalue Problems with 
the Method of Lifting 



S. Setayeshgar ^'*'\ H. B. Keller ^ 

^Department of Applied Mathematics, California Institute of Technology, 
Pasadena, California 91125 

J. E. Pearson^ 

^XCM, MS F645, Los Alamos National Laboratory, Los Alamos, New Mexico 

87545 



Abstract 

Consider an x matrix A for which zero is a defective eigenvalue. In this 
case, the algebraic multiplicity of the zero eigenvalue is greater than the geometric 
multiplicity. We show how an inflated (A^+1) x (A^+1) matrix C can be constructed 
as a rank one perturbation to A , such that C is singular but no longer defective, and 
the nullvectors of C can be easily related to the nullvectors of A. The motivation for 
this construction comes from linear stability analysis of an experimental reaction- 
diffusion system which exhibits the Turing instability. The utility of this scheme 
is accurate numerical computation of nullvector(s) corresponding to a defective 
zero eigenvalue. We show that numerical computations on C yield more accurate 
eigenvectors than direct computation on A. 
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1 Introduction 



In this paper, we present a method for improving the accuracy of computed 
eigenvalues and eigenvectors of defective matrices [1]. This algorithm, which 
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we refer to as lifting, involves construction of a nondefective matrix in a higher 
dimensional space. Standard linear algebra procedures are shown to yield more 
accurate numerical results for the lifted matrix. Our scheme formally resem- 
bles the construction of an augmented system in pseudo-arclength continu- 
ation [2,3,4,5]. Its formulation was motivated by the structure of the linear 
stability analysis of a particular chemical reaction-diffusion system, known as 
the chlorine-dioxide-iodine-malonic-acid system [6,7,8,9] (henceforth referred 
to as CDIMA), which exhibits the Turing instability. 

First, we present the Lifting Theorem. We demonstrate its numerical imple- 
mentation on "small" and "large" defective eigenproblems. We describe the 
connection with the CDIMA linear stability problem, and discuss the general 
applicability of this scheme. 



2 Lifting 

An eigenvalue is defective if its algebraic multiplicity is greater than its ge- 
ometric multiplicity. (The geometric multiplicity of an eigenvalue equals the 
dimension of the space spanned by the eigenvectors corresponding to that 
eigenvalue.) Consider an eigenvalue with geometric multiplicity equal to one: It 
is nondefective if the corresponding left and right eigenvectors are nonorthog- 
onal. For a defective eigenvalue, lifting corresponds to a geometric procedure 
of casting the corresponding eigenvectors of the matrix and the adjoint matrix 
into a higher dimensional space in which these vectors are no longer orthogo- 
nal. In this section, we show how to construct such a well-conditioned matrix 
in the higher dimensional space from a defective matrix. 

For simplicity, we define the notation: 

algebraic multiphcity of eigenvalue fi of matrix M = AM{M; jj) 
geometric multiphcity of eigenvalue ji of matrix M = GM(M; jj). 

Although in the following we consider zero eigenvalues and corresponding 
nuUvectors, the discussion holds generally for nonzero eigenvalues, with A ~ 
M — /il^ , where /at is the appropriate N x N identity matrix. 

Theorem: Consider an x matrix A, such that zero is a defective eigen- 
value of A with AM {A; 0) = n and GM{A; 0) = 1, n > 1. The right and left 
unit nuUvectors of A are given by and i/'"'', respectively. An (A^+1) x (A^+1) 
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matrix C constructed according to 



A 

0^ 



+ vw . 



where 



V 



w 



w 



UJ 



(2) 



will have a simple zero eigenvalue if 

(i) w^(/>^0, 

(ii) V'^v^O, 

(iii) Tjcu 7^ 0. 

V and w are iV-dimensional vectors, and rj and a; are scalars. 
Proof: First, we show there exists a nontrivial vector given by 



(3) 
(4) 
(5) 



(6) 



such that 

= 0. 



(7) 



For $ to be a right nuUvector of C, the N components of {C in the original 
space, as well as the {N + 1)*'* component in the inflated space must be zero: 



A X + av = 

ar] — 0, 

where 

a — w'''x + uj^. 



(8) 
(9) 



(10) 
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Since 77 7^ 0, for Eq. 9 to be satisfied, we must have a = 0. Then, Eq. 8 imphes 
X = r</), where r is a constant, and ^ must be given by 

i = -r—^. 11 



Hence, the pointwise projection of the right nullvector of C onto the origi- 
nal A^-dimensional space is the right nullvector of A, up to a multiplicative 
constant.[3] Additionally, its [N + 1)*'^ component, ^, is nonzero. 

Likewise, we can show that C has a nontrivial left nullvector given by 



* = , (12) 



where 



y = sV, C = -s . (13) 



Again, s is a constant, and we have C 7^ 0- 

Next, we must show that zero is a simple eigenvalue of C This will be done 
in two steps: 



(i) GM(£;0) = 1, 

(ii) = V V + ^ 0. (14) 

First, we prove GM{C; 0) = 1 by contradiction. Assume GM{C; 0) > 1. Then, 
there exists 7^ $ such that 



= 0, 




Eqs. 8 and 9 imply Ax' = 0. We distinguish two cases: 

a. If x' 7^ r(f), then x' is a right nullvector of A, different from r0. This violates 
the assumption GM{A; 0)=1. 

b. If x' = r(f), we must have C,' ^ ^ for 7^ $ . However, from Eqs. 9 and 10, 
^' = —rw^ip/uj = ^. Therefore, is not different from 

^ If $ is normalized to give the unit right nullvector, then = (l + w^^/w) ^ . 
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Hence, wc have shown that GM{C; 0) = GM{A; 0) = 1. (This proof is easily 
generaUzcd to GM{jC;0) = GM{A;0) = v.) To show (ii), we note that zero 
being a defective eigenvalue of A with GM[A; 0) = 1 implies il)^(f) = 0. 
However, the inner product given by Eq. 14 will still be nonzero, since ^ 0. 
This completes the proof that zero is a simple eigenvalue of C 

It is important to consider the generic case in numerical applications, where 
A is not exactly defective but almost-defective, pa 0. For zero to be 

a nondefective eigenvalue of the inner product given by Eq. 14 must be 
nonzero: 



For example, constructing C with v and w such that v = 0, ?7 = — 1, and 
iij — w^(/>, would violate the condition above. Therefore, when A is almost- 
defective, Eq. 16 must be additionally satisfied to guarantee that C will not 
be defective. Below, we give a general prescription for choosing the lifting 
vectors, which unlike the pathological special case just presented, should avoid 
construction of a defective lifted matrix. 

Consider choosing (for simplicity) r] — lo — \^ and (v, w) = (/3vr,7Wr): 

Vrj = rand(i)/iVv, w^^ = rand(j)/iVw, (17) 



where rand(i) are random numbers on [—1, 1], A^v and are normalizing 
factors leading to unit random vectors, and 7 and /? are constants. This con- 
struction should guarantee that when A is singular and exactly defective 
(with AM{A\^) — n, GM{A;0) = 1), or close to defective, then zero will 
be a simple eigenvalue of C Wc will refer to the constants 7 and f3 as lifting 
parameters, since the (A^ -|- 1)*^ components of the left/right nuUvectors of jC 
are proportional to these constants, respectively, and measure the projections 
of these vectors onto the new, infiated direction. 

To extend the above theorem and proof to the case where A is defective and 
AM {A; 0) > GM{A; 0) > 1, we must show that C can be constructed such 
that zero is a nondefective, multiple eigenvalue with AM{C; 0) = GM{C; 0) = 
GM{A; 0). With multiple zero eigenvalues, C will be defective if there exists 
a vector g such that 

^g = *i, (18) 



for some $i G J\f{C), where i = 1 . . .u, and M{C) is the nuUspace of C 
(If such a vector g existed, then it would be a generahzed eigenvector of C, 
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associated with an m x m Jordan block, where m > 1.) Taking the inner 
product of Eq. 18 with *T, where "^J e ^/{C^), j — 1 . . .v, gives 

*J$i = 0. (19) 



Therefore, for C to be nondefective, it must be constructed such that no vector 
in the right nuUspace of C is orthogonal to its entire left nuUspace. 



3 Numerical Results 

3.1 2x2 Defective Eigenproblem 

Consider the following matrix, M: 
( 



M 



TT 1 

\^-7rV4 e 



(20) 



The eigenvalues of M are given by 



7r + e±'^(e2-27re) 



/2. (21) 



The two components of an eigenvector (j)^— {(pi, ^2) satisfy 

02/01 = /X - TT. (22) 



where fi = fi±. For e = 0, M is defective. We use the lifting algorithm to 
explore the error in computing the eigenvector associated with when M 
is close to being defective, e <S 1. The lifted matrix C is constructed with 
A — M — 11^12 according to 



where (vr, Wr) are two-dimensional random unit vectors, and we set 7 = /3 
without loss of generality. 



6 



The right unit nuUvector $ of is computed Tising MATLAB's eigO function 
[10,11]. We explore the hfting error, constructed according to 



£(e,/3) = |$2/$i-(/i+-7r)|. 



(24) 



as a function of parameters e and p. 

Figure 1 shows S{e,j3) as a function of the hfting parameter f3 for different 
values of e. The basic trend is that of decreasing error as a function of increas- 
ing f3 for small to moderate values of /?. In this figure, each point corresponds 
to the mean error computed using 1000 unit random lifting vectors, (vr, Wr). 
The root-mean-square (rms) of the distribution of errors is of the order of the 
mean error for small values of lifting parameter. For /3 >,0(1), the rms is at 
most an order of magnitude greater than the mean error for values of e close 
to machine precision; for larger values of e, the rms is equal to a few times 
the mean error. We have verified that the larger values of rms are a result 
of a few outliers. In Figure 2, for hfting parameter /3 = 1.0, we compare the 
mean error with lifting (open circles) and without lifting (open triangles) as a 
function of e. In this figure, we also show the error using lifting vectors given 
by (v, w)= (?/?, </)), such that the nonorthogonality conditions Eqs. 3 and 4 are 
satisfied by definition (open squares). We note that the lifting error is within 
an order of magnitude of machine precision for all values of e. 

The condition number of the simple singular value of the lifted matrix is given 
by the reciprocal of 



where * and # are the left and right unit nuUvectors oi C In Figure 3, we 
show the condition number as a function of the lifting parameter /3 for different 
values of e, using a single pair of unit random lifting vectors: For e <^ 1, the 
condition number clearly improves with increasing lifting parameter. 

3.2 N X N Defective Eigenproblem 

For applications in which the matrix is a discrete approximation to a continu- 
ous operator, A will have dimension ^ 1. To explore the method of hfting 
in computing the defective eigenvector of larger matrices, we constructed such 
an example from the 2x2 matrix above: with M in the top left block, along 
with the (A^ — 2) x (A^ — 2) block tridiagonal matrix from Poisson's equation 
in the bottom right, we applied a similarity transformation to the resulting 
N X N matrix to obtain a general matrix with no special properties (such 



s(0) = , 



(25) 
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as symmetry or handedness). The eigenvahie. of this matrix is exactly 
known, and the hfted matrix L is constructed as in the 2x2 case. 

Once the right nuUvector of given by = (x, ^) is numerically computed, 
the hfting error is constructed as before: 

£(e,/?) = |(5x)2/(5x)i-(/.+ -7r)|, (26) 

(where S is the similarity transformation). In Figure 4, we show the mean 
error and its rms computed using 50 unit random lifting vectors, (vr, Wr). We 
find that for small values of f3, the lifting error is dominated by the error in 
computing the "zero" eigenvalue, Aq, of jC. For /3 >,C(1), the magnitude of 
this eigenvalue remains at machine precision, while the lifting error is larger. 
For values of /3 such that the magnitude of the "zero" eigenvalue is at ma- 
chine precision, the lifting error is dominated by the error in computing the 
associated eigenvector. We note that there is an "optimal" value of f3 giving 
the smallest lifting error, which is of order unity. For values of (3 much larger 
than unity, the error increases as a function of increasing (3. 



4 CDIMA Eigenproblem 

The symmetry-breaking instability of a system from a homogeneous steady 
state to a patterned state, predicted in 1952 by Alan Turing [12], was observed 
for the first time in the chlorite-iodide-malonic acid CIO2 -1^-MA (CIMA) 
reaction-diffusion system [13,14]. In practice, due to boundary chemical feeds, 
the steady state is not homogeneous, but rather depends on the single space 
coordinate along the feed gradient. A realistic model of the simpler chlorine 
dioxide-iodine-malonic acid CIO2-I2-MA (CDIMA) reaction, which is similar 
to CIMA in terms of its stationary pattern-forming and dynamical behavior, 
was put forth by Lengyel, Rabai, and Epstein (henceforth referred to as LRE) 
[15,16]. In particular, they demonstrated that the fortuitous choice of starch 
as the color indicator for visualizing the patterns in this system provides the 
necessary rescaling of the diffusion coefficient of the activator (iodide) relative 
to the inhibitor (chlorite) species required for this instability to occur. Starch 
is a large molecule that binds to the gel matrix of the reactor and is thus 
effectively immobile. The reversible binding of iodide with starch results in the 
purple starch-triiodide complex and an effective slowing down of the diffusion 
rate of iodide. 

We begin with the general formulation of an [N + l)-component reaction- 
diffusion system with N mobile and a single immobile species (as in the 
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CDIMA system), following [7]: 
dC 



dt 



T>-V'C + J^{U) + C9. 



(27) 



In the evolution equation above, 

/ 





R 

1 



(28) 



U represents the concentrations of mobile species and W that of the immobile 
species, / = f{U; a) describes the reactions among the mobile species with 
control parameters a, and g = g{U,W; (3) describes the reaction between the 
mobile species and the immobile species with control parameters f3. D is the 
diffusion matrix for the mobile subsystem, and is a constant stiochiometry 
vector. We refer the interested reader to previous works for the description of 
the LRE model of the CDIMA system. 

The steady state of this system in the presence of externally imposed boundary 
feed gradients satisfies: 



f{U,;a) + D-^ = 0, (29) 
g{Us,Ws;(3) = (30) 

with boundary conditions Cs(0) = Co and Cs{i) = C^. We wiU focus on the 
physically relevant special case of scalar diffusion, D = D I , where / is 
the identity matrix, as the diffusion constants of ions in solution are compa- 
rable. It is clear from the formulation given by Eq. 27, that the steady state 
of the mobile subsystem is independent of the immobile reactions. However, 
as illustrated in previous works [6,7,8,9], the immobile subsystem affects the 
stability of the mobile system. The stability of the steady state to a mode 
with transverse wavenumber k (k -L z) leads to the following linear operator: 

''f°]-D{k' + ^,)(\°\+Cdg\ (31) 
0^ / ^ M 0^ ' 



where df = fjj is the Jacobian derivative of / with respect to U, and dg = 
{9u,9w) is the Jacobian derivative of g with respect to {U, W), evaluated at 
the steady state. The eigenvalue problem determining the stabihty of Cs{z) to 
the mode with wavenumber k is given by 

A • = Afe*fe (32) 
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where * = (^'^^,'^^\ The space 5" on which the hnear operator acts can 
be written as S = Sm © 5'/, where Sm and Si are mobile and immobile 
subspaces, respectively. In gradient systems, the basis set spanning Sm is given 
by the Cartesian product of orthogonal unit vectors {ej}, i = 1,. . . ,N, 
corresponding to the N mobile species, and basis functions, defined on the 
interval in the -^-direction, for example given by {sin rmrz/Lz} for Dirichlet 
boundary conditions, where m = 1, . . . ,Nz. Similarly, Sj would be spanned 
by eAT+i^jsinrnTT^/L^}. We make the physical assumption that no instability 
occurs in the immobile subsystem, ^'^^ 7^ 0. 

A Turing bifurcation occurs when a single goes through zero and becomes 
positive as the parameter ot goes through the critical value ac- These condi- 
tions 



dXk 



dk 



(33) 



are satisfied at the onset of the instability, and define the critical wavenumber 
kc- Note that in the case of scalar diffusion, no such instability can occur in 
the absence of reaction with the immobile species. Prom the above. 



0, 



(34) 



and 



dC 

'dk 



'dk 



(35) 



Taking the inner product of Eq. 35 on the left by ^l^, the nuUvector of Cl^, 
we find 



(36) 



for kc 7^ 0. In the absence of the immobile reaction, this implies that zero is a 
defective eigenvalue of C^^ at the Turing bifurcation. However, it follows that 



(37) 



since there must exist a nontrivial solution to 
gv^Z + 9wK = 0- 



(3J 
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(Otherwise a nontrivial null vector would not exist.) Hence, zero is not a 
defective eigenvalue of Ck^- Furthermore, it can be easily shown that the sin- 
gularity of Ck as a function of parameters ot and wavenumber k is independent 
of the immobile reaction [7] . 

In this way, the structure of the linear stability of the {N + l)-component 
reaction-diffusion system at the Turing bifurcation motivated the formulation 
of a general scheme for appropriately lifting a defective eigenvalue problem 
to a higher dimensional space in which it is no longer defective. Despite this 
connection, our numerical examples of the lifting algorithm in the previous 
section did not include linear stability analysis of the CDIMA system. This is 
because the improvement in accurate calculation of the eigenvector associated 
with the defective eigenvalue is best illustrated when the eigenvalue is exactly 
known. For the gradient CDIMA system, the defective eigenvalue jic = Dk^ 
is not known apriori. 



5 Concluding Renictrks 

The stability analysis of the experimental CDIMA reaction-diffusion system 
has led to the formulation of a new scheme in computational linear algebra for 
more accurate numerical solution of defective eigenvectors. We have demon- 
strated improvement in numerical accuracy by several orders of magnitude 
using illustrative toy examples. This scheme should be of utility in physical 
applications where accurate eigenvectors of defective matrices are sought, as 
well as of general use in improving accuracy of matrix computations involving 
such matrices. 
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Fig. 1. The lifting error, £{€,(3), is plotted as a function of lifting parameter, (3, for 
different values of e. The lifting vectors are given by: v = (/3 Vr, 1) and w = (/? Wr, 1), 
where Vr and Wr are two-dimensional unit random vectors. Each point corresponds 
to the mean error obtained using 1000 unit random lifting vectors, (vr, Wr). 
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Fig. 2. Error as a function of e with lifting parameter (3 = 1.0: Open circles show the 
error using lifting vectors given by v = (/3 Vr, 1) and w = (/3wr, 1), where Vr and 
Wr are two-dimensional unit random vectors. Each point corresponds to the mean 
error obtained using 1000 unit random lifting vectors, (vr, Wr). Open squares show 
the error using lifting vectors given by v = {(3tj}, 1) and w = {13(f), 1), where ■0 and 
(f) are the left and right unit nullvectors of A, respectively. The error without lifting 
is given by the long-dashed line (open triangles). 
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Fig. 3. The condition number for the singular value of the 3x3 lifted matrix, C , is 
plotted as a function of hfting parameter, /3, for different values of e. A single, unit 
random lifting vector, (vr, Wr), was used to construct the lifted matrix. 
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Fig. 4. The lifting error (open circles) and the magnitude of the "zero" eigen- 
value, |Ao|, of C (filled squares) axe plotted as a function of lifting parameter, 
/3, with e = 10-12 ^nd N = 500. The error without lifting is shown for comparison 
(long-dashed line). The lifting vectors are given by: v = (/3 Vr, 1) and w = (/3 Wr, 1), 
where Vr and Wr are {N — l)-dimensional unit random vectors. Each point cor- 
responds to the average quantity obtained using 50 unit random lifting vectors, 
(vr, Wr), and the error bars give the rms of the distribution of errors. 
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